function [stored_hist stored_diams labeledgloms] = DepthMeasures_v2(depth,ThresholdRange,DiameterRangePixelUnits)
% function to calculate threshold statistics for glomeruli blobs in a volume image.
%
%  This second version of DepthMeasures uses bwconncomp rather than bwlabeln as the
%  primary connected component algorithm. The binary algorithm uses different slices
%  through the 'depth' volume to collect connected regions.The connected component
%  labels are organised and counted.
%  Input arguments:
%   depth : The input volume. The volume is normalised automatically then
%   used to test for connected components at different threshold levels. 
%   ThresholdRange (optional):    linearly space vector between 0 and 1 that decides which  
%      threshold values are used in the threshold statistics.
%   DiameterRangePixelUnits (optional):  linearly spaced vector controlling the estimated 
%      diameter of labelled components (in pixel units).
%  Outputs:
%   stored_hist: array of histograms at each threshold showing number of voxels per labelled components
%   stored_diams: array of histograms for estimated diameters of the labelled components
%   labelledgloms: volume thresholded at midpoint with labelled glomeruli 
%
%   - Michael Eager,  (michael.eager@monash.edu)
%   - (c) 2012,  Monash Biomedical Imaging,  Monash University, Australia

%     Copyright © 2012-2013 Michael Eager <michael.eager@monash.edu> 
%
%     This file is part of Xglom.
% 
%     This is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.


display('Normalising depth image.');
depth_min=min(depth(:));
depth_max=max(depth(:));
if depth_min ~= 0 && dept_max ~= 1
    depth = depth + abs(depth_min);
    depth = depth./max_depth;
end

if nargin < 2
ThresholdRange= 0.01:0.01:0.5;
end

if nargin < 3
DiameterRangePixelUnits = 1:0.1:15;
end

stored_hist=zeros(50,1000);
stored_diams = zeros(50,numel(DiameterRangePixelUnits));
counter=10;
%% Thresholding and discard single voxels
tic
for counter=1:length(ThresholdRange) % Loop over different Thresholds
    %%
    %Thresholding
    Holes=depth>ThresholdRange(counter);
    
    %Discard single voxels
    %     ConArray=conndef(3,'minimal');
    %     ConArray(2,2,2)=0;
    %     Holes=imdilate(Holes,ConArray).*Holes;
    %     Holes=single(Holes);
    
    %% Statistics
    %Label Glomerular Regions
    tic;CC = bwconncomp(Holes,6);toc
    
    if counter == 30
        labeledgloms = labelmatrix(CC);
    end
    
    %calculate area histogram
    areas = cellfun(@numel,CC.PixelIdxList);
    area_max = 1000; %max(areas);
    area_range=(1:area_max);
    [area_hist, area_axis]=hist(areas,area_max);
    stored_hist(counter,:) = area_hist;
    
    %Total number of detected objects in area_range
    ['sum(area_hist(12:100)):  ', num2str(sum(area_hist))];
    
    %Diameter histogram from area data assuming spherical glomeruli
    [h x]=hist((6/pi*areas).^(1/3),DiameterRangePixelUnits); h(end)=0;
    stored_diams(counter,:) = h;

end
toc




%% Find Centroids of labels
% locations = regionprops(D,'centroid');
% glom_centroids = cat(1,locations.Centroid);
% xloc = glom_centroids(:,1);
% yloc = glom_centroids(:,2);
% zloc  = glom_centroids(:,3)

% OR
%% STATS = REGIONPROPS(CC,PROPERTIES);  % where CC is the ans of bwconncomp

%     Shape Measurements
%
%       'Area'              'EulerNumber'       'Orientation'
%       'BoundingBox'       'Extent'            'Perimeter'
%       'Centroid'          'Extrema'           'PixelIdxList'
%       'ConvexArea'        'FilledArea'        'PixelList'
%       'ConvexHull'        'FilledImage'       'Solidity'
%       'ConvexImage'       'Image'             'SubarrayIdx'
%       'Eccentricity'      'MajorAxisLength'
%       'EquivDiameter'     'MinorAxisLength'
%
%     Pixel Value Measurements (requires grayscale image as an input)
%
%       'MaxIntensity'
%       'MeanIntensity'
%       'MinIntensity'
%       'PixelValues'
%       'WeightedCentroid'
%



